Когда признаков много, можно представить все объекты как облако точек в многомерном пространстве.
Migration by Don McCullough on Flickr
black shadows for a white horses / les negres ombres dels cavalls blancs by Ferran Jordà on Flickr
Nonmetric Multidimensional Scaling, nMDS (Shepard 1962, 1966, Kruskal 1964a, b)
Метод визуализации отношений между объектами в пространстве с небольшим числом измерений (обычно 2).
Исходные данные — матрица расстояний между объектами в многомерном пространстве.
nMDS ординация старается сохранить отношения между объектами.
Расстояния по автодорогам не совсем “евклидовы”: например, из Кронштадта до Санкт-Петербурга нельзя добраться по прямой.
Карта
Проекция при помощи анализа главных координат (PCoA)
Карта
Ординация nMDS
Ординация сохраняет ранг расстояний между объектами (похожие располагаются близко, непохожие — далеко; если объект А похож на B, больше чем на C, то и на ординации он, скорее всего, окажется ближе к B, чем к C).
На ординации nMDS имеет смысл только взаиморасположение объектов. Облако точек в осях MDS можно вращать, перемещать, зеркально отражать. Суть ординации от этого не изменится.
Численные значения координат точек в ординации лишены смысла (их вообще можно не приводить на итоговой ординации).
Взаиморасположение точек на плоскости подобно взаиморасположению точек в многомерном пространстве признаков, но значения не полностью совпадают.
Диаграмма Шепарда — график зависимости расстояний между точками на ординации и соответствующих им значений коэффициентов различия в исходном пространстве.
В идеале, ранги расстояний между точками \(d_{h,i}\) на хорошей ординации nMDS должны совпадать с рангами коэффициентов различия \(D_{h,i}\) в исходном многомерном пространстве.
Форма и расположение облака точек на диаграмме Шепарда позволяет судить о качестве ординации.
В идеале, большим расстояниям в исходном пространстве должны соответствовать большие расстояния на ординации, а меньшим — меньшие (т.е. должны сохраняться ранги расстояний).
Красная линия — предсказанные монотонной регрессией значения расстояний на ординации — монотонно возрастает.
Монотонная регрессия минимизирует сумму квадратов отклонений значений координат на ординации от значений в исходном многомерном пространстве, сохраняя монотонные отношения между ними.
Алгоритм подбора монотонной регрессии называется Pool-Adjacent-Violator-Algorithm (PAVA)
Стресс (Stress) показывает, насколько соответствуют друг другу взаиморасположение точек на плоскости ординации и в исходном многомерном пространстве признаков.
Стресс 1 (Stress(1), Kruskal, 1964a) — это корень из суммы квадратов остатков монотонной регрессии, отнесенной к сумме квадратов всех коэффициентов различия в исходной матрице.
\[ Stress(1) = \sqrt{\frac{\sum_{h < i}{(d_{h,i} - \hat d_{h,i})^2}}{\sum_{h < i}{d_{{h,i}}^2}}}\]
Хорошо, когда \(Stress(1) < 0.2\).
Данные можно скачать с сайта Pangaea (Krapivin 2017)
dat <- read.delim("data/Krapivin_2017_Medulis-symb.tab", skip = 36,
stringsAsFactors = TRUE)
str(dat)
## 'data.frame': 548 obs. of 23 variables: ## $ Event : Factor w/ 3 levels "Gorelyi_Isl",..: 3 3 3 3 3 3 3 3 3 3 ... ## $ Date.Time : Factor w/ 1 level "2013-07": 1 1 1 1 1 1 1 1 1 1 ... ## $ Ord.No : int 1 2 3 4 5 6 7 8 9 10 ... ## $ Area : Factor w/ 3 levels "B Solovkecki Isl",..: 1 1 1 1 1 1 1 1 1 1 ... ## $ Site : int 1 1 1 1 1 1 1 1 1 1 ... ## $ Zone : Factor w/ 3 levels "down","middle",..: 1 1 1 1 1 1 1 1 1 1 ... ## $ Comment : int 1 1 1 1 1 1 1 1 1 1 ... ## $ M..edulis.L.shell..mm.: num 31.2 41.1 46.4 40 34.5 47.6 34.3 49.6 49.8 30 ... ## $ M..edulis.age..a. : num 1 2 2 3 4 3 2 4 3 1 ... ## $ Urastoma.spp..... : int 26 16 19 11 14 16 10 4 11 0 ... ## $ Renicola.spp..... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Himasthla.spp..... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Metacercaria.... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Gymnophallus.spp..... : int 0 0 0 0 0 0 0 0 4 0 ... ## $ Bucephalidae.... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Alg.... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Nematoda.... : int 1 1 0 0 0 0 0 0 0 0 ... ## $ Microsetella.... : int 0 0 0 0 0 2 0 0 0 0 ... ## $ Copepoda.... : int 0 0 0 0 0 1 0 0 0 0 ... ## $ Chironomidae.... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Halacaridae.... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Jaera.spp..... : int 0 0 0 0 0 0 0 0 0 0 ... ## $ Ostrac.... : int 0 0 0 0 0 0 0 0 0 0 ...
colnames(dat) # Было
## [1] "Event" "Date.Time" "Ord.No" ## [4] "Area" "Site" "Zone" ## [7] "Comment" "M..edulis.L.shell..mm." "M..edulis.age..a." ## [10] "Urastoma.spp....." "Renicola.spp....." "Himasthla.spp....." ## [13] "Metacercaria...." "Gymnophallus.spp....." "Bucephalidae...." ## [16] "Alg...." "Nematoda...." "Microsetella...." ## [19] "Copepoda...." "Chironomidae...." "Halacaridae...." ## [22] "Jaera.spp....." "Ostrac...."
colnames(dat) <- gsub("[.]", replacement = "", colnames(dat))
dat <- rename(dat, L = "MedulisLshellmm", Age = "Medulisagea")
colnames(dat) # Стало
## [1] "Event" "DateTime" "OrdNo" "Area" ## [5] "Site" "Zone" "Comment" "L" ## [9] "Age" "Urastomaspp" "Renicolaspp" "Himasthlaspp" ## [13] "Metacercaria" "Gymnophallusspp" "Bucephalidae" "Alg" ## [17] "Nematoda" "Microsetella" "Copepoda" "Chironomidae" ## [21] "Halacaridae" "Jaeraspp" "Ostrac"
# Делаем сайт фактором
dat$Site <- factor(dat$Site, levels = c(2, 3, 1), labels = c("G","L","S"))
# Сливаем редкие виды в одну категорию
f_remove <- c("Nematoda", "Microsetella", "Copepoda", "Chironomidae",
"Halacaridae", "Jaeraspp", "Ostrac")
dat$Other <- rowSums(dat[, f_remove])
# Суммарная численность симбионтов
f_sp <- c("Urastomaspp", "Renicolaspp", "Himasthlaspp", "Metacercaria",
"Gymnophallusspp", "Alg", "Other")
dat$Total <- rowSums(dat[, f_sp])
# Данные для анализа
# Только мидии с симбионтами и возрастом от 3 до 8 лет
dat <- dat[dat$Total != 0 & dat$Age %in% 3:8, ]
spec <- dat[, f_sp] # виды-симбионты
env <- dat[, c("Zone", "Site", "L", "Age")] # свойства мидий-хозяев
Постройте ординацию мидий по обилию разных видов-симбионтов с использованием коэффициента различия Брея-Куртиса. Следите, чтобы алгоритму удалось найти финальное решение. Если необходимо, настройте вызов metaMDS()
Вычислите стресс для получившейся ординации.
Нарисуйте простейший график.
ord_mus <-
ord_mus <- metaMDS(spec, dist = "bray")
По-умолчанию metaMDS() пытается подобрать оптимальную трансформацию, подходящую для анализа сообществ. В данном случае — извлекает корень и делает двойную висконсинскую стандартизацию.
Square root transformation Wisconsin double standardization Run 0 stress 0.1319177
Такая модель не сходится.
*** No convergence -- monoMDS stopping criteria:
12: stress ratio > sratmax
8: scale factor of the gradient < sfgrmin
После отключения автоматического подбора трансформации данных удается найти решение.
ord_mus <- metaMDS(spec, dist = "bray", autotransform = FALSE)
## Run 0 stress 0.141 ## Run 1 stress 0.145 ## Run 2 stress 0.144 ## Run 3 stress 0.147 ## Run 4 stress 0.142 ## Run 5 stress 0.141 ## ... New best solution ## ... Procrustes: rmse 0.00132 max resid 0.0123 ## Run 6 stress 0.163 ## Run 7 stress 0.153 ## Run 8 stress 0.141 ## ... Procrustes: rmse 0.00181 max resid 0.0326 ## Run 9 stress 0.141 ## ... Procrustes: rmse 0.000133 max resid 0.002 ## ... Similar to previous best ## Run 10 stress 0.141 ## ... Procrustes: rmse 0.00339 max resid 0.0334 ## Run 11 stress 0.155 ## Run 12 stress 0.141 ## ... Procrustes: rmse 0.00121 max resid 0.0119 ## Run 13 stress 0.141 ## ... Procrustes: rmse 0.0014 max resid 0.0119 ## Run 14 stress 0.156 ## Run 15 stress 0.156 ## Run 16 stress 0.141 ## ... New best solution ## ... Procrustes: rmse 0.00062 max resid 0.0093 ## ... Similar to previous best ## Run 17 stress 0.141 ## ... Procrustes: rmse 0.00365 max resid 0.0337 ## Run 18 stress 0.146 ## Run 19 stress 0.141 ## ... Procrustes: rmse 0.00235 max resid 0.0335 ## Run 20 stress 0.146 ## *** Solution reached
Хорошее решение, но график нужно украсить.
ord_mus$stress
## [1] 0.141
ordiplot(ord_mus)
# Палитры
pal_col <- c("red", "green", "steelblue")
pal_sh <- c(1, 2, 0)
# Украшенный график
ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])
ordiplot(ord_mus, type = "n") points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site]) text(ord_mus, display = "spec", cex = 0.9, col = "grey20")
Используя данные, приведенные ниже, постройте график ординации при помощи ggplot2. Покажите точками — мидий, текстом — центроиды для видов-симбионтов, цветом — зону литорали, формой — сайт, размером маркеров — размер мидий.
# Координаты точек (мидий) ord_mus_pt <- data.frame(env, scores(ord_mus, display = "sites")) head(ord_mus_pt, 2)
# Координаты центроидов переменных (видов-симбионтов) ord_mus_sp <- data.frame(scores(ord_mus, display = "species")) ord_mus_sp$Species <- rownames(ord_mus_sp) head(ord_mus_sp, 2)
gg_ord_mus <- ggplot() +
geom_point(data = ord_mus_pt,
aes(x = NMDS1, y = NMDS2, colour = Zone,
shape = Site, size = L), alpha = 0.5)
gg_ord_mus
gg_ord_mus_sp <- gg_ord_mus +
geom_text(data = ord_mus_sp,
aes(x = NMDS1, y = NMDS2,
label = Species))
gg_ord_mus_sp
envfit()Функция envfit() строит регрессионную модель зависимости переменной среды от координат точек на плоскости ординации.
Модель вот такого вида:
\(y_i = \beta_1 \text{NMDS1}_i + \beta_2 \text{NMDS2}_i + \varepsilon_i\)
Значимость зависимости оценивается при помощи пермутационного теста (это позволяет работать в ситуации, когда не полностью выполнены условия применимости).
Такой вызов функции позволит подобрать сразу несколько моделей и изобразить их на графике ординации:
ef <- envfit(ord_mus, env[, c("Zone", "Site", "L", "Age")])
envfit() для непрерывных переменныхef$vectors
## NMDS1 NMDS2 r2 Pr(>r) ## L -1.000 -0.022 0.19 0.001 *** ## Age 0.777 -0.629 0.04 0.001 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## Permutation: free ## Number of permutations: 999
MDS1 и MDS2 содержат информацию о направлении вектора непрерывной переменной. Длина вектора пропорциональна корню из коэффициента детерминации для модели.r2 — \(R^2\) из регресcионного анализа.Pr(>r) — пермутационная оценка статистической значимости.envfit() для дискретных переменныхef$factors
## Centroids: ## NMDS1 NMDS2 ## Zonedown -0.65 -0.09 ## Zonemiddle -0.14 -0.05 ## Zoneup 0.83 0.14 ## SiteG 0.19 0.17 ## SiteL 0.23 -0.41 ## SiteS -0.42 0.23 ## ## Goodness of fit: ## r2 Pr(>r) ## Zone 0.37 0.001 *** ## Site 0.16 0.001 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## Permutation: free ## Number of permutations: 999
MDS1 и MDS2 содержат координаты центроидов для дискретных факторов (т.е. положение центра облака точек, относящихся к категории фактора).r2 — \(R^2\) из регресcионного анализа.Pr(>r) — пермутационная оценка статистической значимости.envfit()ordiplot(ord_mus, type = "n") points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site]) plot(ef)
# install.packages("devtools")
# devtools::install_github("gavinsimpson/ggvegan")
library(ggvegan)
ord_mus_ef <- fortify(ef)
ord_mus_ef
envfit()gg_ord_mus +
geom_segment(data = ord_mus_ef[ord_mus_ef$Type == "Vector", ],
aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
arrow = arrow(length = unit(0.25, "cm"))) +
geom_text(data = ord_mus_ef[ord_mus_ef$Type == "Vector", ],
aes(x = NMDS1, y = NMDS2, label = Label, hjust = 1.1, vjust = 1)) +
geom_text(data = ord_mus_ef[ord_mus_ef$Type == "Centroid", ],
aes(x = NMDS1, y = NMDS2, label = Label, hjust = 1.1, vjust = 1))
envfit()envfit() основан на предположении о линейном характере связи между внешней переменной и координатами на ординации. Это не всегда так.ordisurf()ordisurf() применяется, когда зависимость нелинейна.
В основе работы функции — общая аддитивная модель (General Additive Model, GAM).
GAM позволяют подобрать нелинейные зависимости, даже если заранее неизвестна их форма.
ordisurf() использует пакет mgcv для подбора GAM на основе метода тонких пластин (thin plate spline).
ordisurf()Изолинии на графике ординации показывают предсказанное значение зависимой переменной для разных участков ординации.
par(mfrow = c(1, 2)) os_L <- ordisurf(ord_mus, env$L, method = "REML") os_Age <- ordisurf(ord_mus, env$Age, method = "REML") par(mfrow = c(1, 1))
ordisurf()Длина мидий линейно меняется на ординации (эффективное число степеней свободы для сплайна \(edf=1.96\)). Зависимость довольно хорошо выражена (объясненная девианса 18.7%).
summary(os_L)
## ## Family: gaussian ## Link function: identity ## ## Formula: ## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE) ## ## Parametric coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 43.410 0.484 89.7 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: ## edf Ref.df F p-value ## s(x1,x2) 1.96 9 9.83 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## R-sq.(adj) = 0.183 Deviance explained = 18.7% ## -REML = 1457.6 Scale est. = 92.442 n = 395
ordisurf()Возраст мидий меняется нелинейно вдоль ординационных осей (эффективное число степеней свободы для сплайна \(edf=4.99\)). Эта зависимость довольно слабая (объясненная девианса 8.97%).
summary(os_Age)
## ## Family: gaussian ## Link function: identity ## ## Formula: ## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE) ## ## Parametric coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 4.8861 0.0714 68.5 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Approximate significance of smooth terms: ## edf Ref.df F p-value ## s(x1,x2) 5.01 9 3.73 0.00000097 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## R-sq.(adj) = 0.0785 Deviance explained = 9.02% ## -REML = 704.9 Scale est. = 2.0111 n = 395
fortify_ordisurf <- function(model) {
# Fortifies an object of class ordisurf to produce
# a dataframe of contour coordinates
# for subsequent plotting in ggplot.
xy <- expand.grid(x = model$grid$x, y = model$grid$y)
xyz <- cbind(xy, c(model$grid$z))
names(xyz) <- c("x", "y", "z")
return(na.omit(xyz))
}
ord_mus_os <- fortify_ordisurf(os_Age)
head(ord_mus_os, 4)
ordisurf()ggplot(data = ord_mus_os, aes(x = x, y = y, z = z)) + stat_contour(aes(colour = ..level..)) + labs(x = "NMDS1", y = "NMDS2", colour = "Age")
ordisurf()ordisurf() важно, чтобы в ординации участвовало много объектов и они располагались на плоскости ординации более-менее равномерно. Пустоты и отскоки снижают качество предсказаний.
envfit() и ordisurf() обращайте внимание на форму связи.Во всех примерах необходимо:
Источники данных
mafragh, пакет ade4).BCI, пакет vegan).